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ABSTRACT 

The long-term goal initiated in this work is to obtain fast 
algorithms and implementations for definite integration in 
Almkvist and Zeilberger's framework of (differential) cre- 
ative telescoping. Our complexity-driven approach is to ob- 
tain tight degree bounds on the various expressions involved 
in the method. To make the problem more tractable, we re- 
strict to bivariate rational functions. By considering this 
constrained class of inputs, we are able to blend the general 
method of creative telescoping with the well-known Hermite 
reduction. We then use our new method to compute diago- 
nals of rational power series arising from combinatorics. 

Categories and Subject Descriptors: 
1.1.2 [Computing Methodologies]: Symbolic and Alge- 
braic Manipulations — Algebraic Algorithms 

General Terms: Algorithms, Theory. 

Keywords: Hermite reduction, creative telescoping. 

1. INTRODUCTION 

The long-term goal of the research initiated in the present 
work is to obtain fast algorithms and implementations for 
the definite integration of general special functions, in a 
complexity-driven perspective. 

As most special-function integrals cannot be expressed in 
closed form, their evaluation cannot be based on table look- 
ups only, and even when closed forms are available, they 
may prove to be intractable in further manipulations. In 
both cases, the difficulty can be mitigated by representing 
functions by annihilating differential operators. This moti- 
vated Zeilberger to introduce a method now known as cre- 
ative telescoping [20], which applies to a large class of special 
functions: the D-finite functions [14] defined by sets of linear 
differential equations of any order, with polynomial coeffi- 
cients. Zeilberger's method applies in general to multiple 
integrals and sums. 
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A sketch of Zeilberger's method is as follows. Given a D- 
finite function / of the variables x and y, the definite integral 

F( x ) ~ Sa f( x > y) dy i s D-finite, and a linear differential 
equation satisfied by F can be constructed [20]. To explain 
this, let k be a field of characteristic zero, D x and D v be 
the usual derivations on the rational-function field k(x,y), 
both restricting to zero on k, and let k(x, y){D x , D y ) be the 
ring of linear differential operators over k(x, y). The heart of 
the method is to solve the differential telescoping equation 
(1) below for L £ k[x](D x ) \ {0} and g = R(f) for some 
R £ k(x,y){D x , D y ) . The operator L is called a telescoper 
for /, and g a certificate of L for /. Under the assumption 

lim g(x,y) — lim g(x,y) for x in some domain, 

L(x, D x ) is then proved to be an annihilator of F. 

The main emphasis in works since the 1990's has been 
on finding telescopers of order minimal over all telescopers 
for /, which are called minimal telescopers. (Two minimal 
telescopers differ by a multiplicative factor in k(x).) In view 
of the computational difficulty of solving (1), there has been 
special attention to subclasses of inputs. Of particular im- 
portance is the case of hyperexponential functions, defined 
by first-order differential equations, studied by Almkvist and 
Zeilberger in [1]. Their method is a direct differential ana- 
logue of Zeilberger's algorithm for the recurrence case [21]. 

On the other hand, very little is known about the com- 
plexity of creative telescoping: the only related result seems 
to be an analysis in [9] of an algorithm for hyperexponen- 
tial indefinite integration. In order to get complexity esti- 
mates, we simplify the problem by restricting to a smaller 
class of inputs, namely that of bivariate rational functions. 
Although restricted, this class already has many applica- 
tions, for instance in combinatorics, where many nontrivial 
problems are encoded as diagonals of rational formal power 
series, themselves expressible as integrals. Our goal thus 
reads as follows. 

Problem Given f = P/Q € k(x, y) \{0}, find a pair (L, g) 
with L — ^2i^ r/i(x)D x in k[x]{D x ) \ {0} and g in k(x,y) 
such that 

L(x,D x )(f) = D y (g). (1) 

By considering this more constrained class of inputs, we 
are indeed able to blend the general method of creative tele- 
scoping with the well-known Hermite reduction [10]. 

Essentially two algorithms for minimal telescopers can be 
found in the literature: The classical way [1] is to apply a 
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Figure 1: Complexity of creative telescoping methods (under Hyp. (H')), together with bounds on output 



differential analogue of Gosper's indefinite summation algo- 
rithm, which reduces the problem to solving an auxiliary 
linear differential equation for polynomial solutions. An al- 
gorithm developed later in [7] (see also [12]) performs Her- 
mite reduction on / to get an additive decomposition of the 
form / = D y (a) + Xl"=i u i/ v i, where the it; and Vi are in 
k(x)[y] and the Vi are squarefree. Then, the algorithm in [1] 
is applied to each Ui/vt to get a telescoper Li minimal for it. 
The least common left multiple of the Li's is then proved to 
be a minimal telescoper for /. This algorithm performs well 
only for specific inputs (both in practice and from the com- 
plexity viewpoint), but it inspired our Lemma 22 via [12]. 

As a first contribution in this article, we present a new, 
provably faster algorithm for computing minimal telescopers 
for bivariate rational functions. Instead of a single use of 
Hermite reduction as in [12] , we apply Hermite reduction to 
the D x (f)'s, iteratively for i — 0,1, ... , which yields 

D i x (f) = D y (g i ) + ^ (2) 

for some factor w of the squarefree part of the denominator 
of /. If 770, . • . , Tj p £ k(x) are not all zero and such that 
Yli=o ViWi ~ 0, then the operator X]f=o Vi^x i s a telescoper 
for /, and more specifically, the first nontrivial linear relation 
obtained in this way yields a minimal telescoper for /. 

As a second contribution, we give the first proof of a poly- 
nomial complexity for creative telescoping on a specific class 
of inputs, namely on bivariate rational functions. For min- 
imal telescopers, only a polynomial bound on d x (but none 
on d y ) was given for special inputs in [7]; more specifically, 
we derive complexity estimates for all mentioned methods 
(see Fig. 1), showing that our approach is faster. Further- 
more, we analyse the bidegrees of non minimal telescop- 
ers generated by other approaches: Lipshitz' work [13] can 
be rephrased into an existence theorem for telescopers with 
polynomial size; the approach followed in the recent work on 
algebraic functions [3] leads to telescopers of smaller degree 
sizes. These are new instances of the philosophy, promoted 
in [3], that relaxing minimality can produce smaller outputs. 

A third contribution is a fast Maple implementation [22], 
incorporating a careful implementation of the original Her- 
mite reduction algorithm, making use of the special form 
oiwi/w in (2) and of usual modular techniques (probabilis- 
tic rank estimate) to determine when to invoke the solver 
for linear algebraic equations. Experimental results indicate 
that our implementation outperforms Maple's core routine. 

Note that for the fastest method we propose, denoted 
by HI in Tables 1-3, we chose to output the certificate as 
a mere sum of (small) rational functions, without any form 
of normalisation. This choice seems to be uncommon for 
creative-telescoping algorithms, but a motivation is how the 
certificate is used in practice: Very often, like for appli- 
cations to diagonals in § 5, the certificate is actually not 
needed. In other applications, the next step of the method 



of creative telescoping is to integrate (1) between a and (3, 
leading to L(F)(x) — g(x, a) — g(x, (3). Therefore, only eval- 
uations of the certificate are really needed, and normalisa- 
tion can be postponed to after specialising at 01 and /3. 

The end of this section, § 1.1, provides classical complexity 
results, notation, and hypotheses that will be used through- 
out. We then study Hermite reduction over k(x) in §2, 
proving output degree bounds and a low-complexity algo- 
rithm. This is then applied in § 3 to derive our new algo- 
rithm for creative telescoping, and to compare its complexity 
with that of Almkvist and Zeilberger's approach. For non- 
minimal telescopers, we show the existence of some of lower 
arithmetic size in § 4: cubic for nonminimal order instead 
of quartic for minimal order. See the summary in Figure 1, 
where the low complexity of algorithms for minimal tele- 
scopers relies on Storjohann and Villard's algorithms [19], 
thus inducing a certified probabilistic feature. We apply our 
results to the calculation of diagonals in § 5, and describe our 
implementation and comment on execution timings in § 5. 

1.1 Background on complexity — Notation 

We recall basic notation and complexity facts for later use. 
Let k be again a field of characteristic zero. Unless other- 
wise specified, all complexity estimates are given in terms of 
arithmetical operations in k, which we denote by "ops". Let 
k[x]™j n be the set ofmxn matrices with coefficients in k[x] 
of degree at most d. Let uj € [2, 3] be a feasible exponent of 
matrix multiplication, so that two matrices from k nxn can 
be multiplied using C(n") ops. Facts 1 and 2 below show the 
complexity of multipoint evaluation, rational interpolation, 
and algebraic operations on polynomial matrices using fast 
arithmetic, where the notation O(-) indicates cost estimates 
with hidden logarithmic factors [6, Def. 25.8]. 

Fact 1 For p £ k[x] of degree less than n, pairwise distinct 
wo, . . . , u n -i in k, and vo, . . . , f n -i £ k, we have: 

(i) Evaluating p at the Ui 's takes 6(n) ops. 
(ii) For m £ {1, . . . , n}, constructing f = s/t £ k(x) with 
deg a ,(s) < m and Aeg x (t) < n — m such that t(ui) 7^ 
and f(ui) = Vi for < i < n — 1 takes 0(n) ops. 

Fact 2 For M in k[x\™* n , d > 0, we have: 

(i) If M = (Mi M2) is an mvertible n x n matrix with 
Mi € k[x]^ j ni , where i = 1, 2 and n\ + n2 = n, then 
the degree of det(M) is at most nidi + n 2 d 2 . 

(ii) If M — (Mi M 2 ) is not of full rank and with Mi £ 
k[x]™£ n ' , where i = 1, 2 and ni + n 2 = n, then there 
exists a nonzero u £ k[x] n with coefficients of degree 
at most nidi + n 2 d 2 such that Mu — 0. 

(in) The rank r and a basis of the null space of M can be 
computed using 0(nmr u '~ 2 d) ops. 



(For proofs, see [6, Cor. 10.8, 5.18, 11.6] and [19, Th. 7.3].) 

We call squarefree factorisation of Q G k[x, y]\k[x] w.r.t. y 
the unique product qQiQ\ ■ ■ ■ Q™ equal to Q for q G k[x] and 
Qi G k[x,y] satisfying deg y (Qm) > and such that the Qi's 
are primitive, squarefree, and pairwise coprime. The square- 
free part Q* of Q w.r.t. y is the product Q1Q2 ■ ■ ■ Qm- Let 
Q~ denote the polynomial Q/Q*, and lc y (Q) the leading 
coefficient of Q w.r.t. y. The following two formulas about 
Q, Q* , and Q~ can be proved by mere calculations. 

Fact 3 Let Qi denote Q* /Qi. Then we have 

(I) Q*D y {Q-)/Q~ = YZxtt ~ l)QiD y (Qi) G k[x,y]; 

(II) D y (Q)/Q~ = J2iLi iQiD y {Qi) G k[x,y]. 

Let / = P/Q be a nonzero element in k(x, y), where P, Q 
are two coprime polynomials in k[x,y]. The degree of / 
in x is defined to be max{deg a .(P), deg x (Q)}, and denoted 
by deg x (f). The degree of / in y is defined similarly. The 
bidegree of / is the pair (deg x (/), deg (/)), which is denoted 
by bideg(/). The bidegree of / is said to be bounded (above) 
by (ct,/3), written bideg(/) < (ct,/3), when deg x (f) < a 
and deg„(/) < p. 

We say that / = P/Q is proper if the degree of P in y 
is less than that of Q. For creative telescoping, we may 
always assume w.l.o.g. that / = P/Q is proper. If not, 
rewrite / = D y (p) + / with p G fc(a;)[j/] and / proper. A 
telescoper L for / with certificate g is a telescoper for / with 
certificate L(p) + g. 

Hypothesis (H) From now on, P and Q are assumed to be 
nonzero polynomials in k[x, y] such that deg y (P) < deg y (Q), 
gcd(P,Q) — 1, and Q is primitive w.r.t. y. 

Notation From now on, we write (d^jdy), (d x ,d*), and 
(d~,d~) for the bidegrees of Q, Q* , and Q , respectively. 

The following hypothesis makes our estimates concise. 

Hypothesis (H') Occasionally, we shall require the ex- 
tended hypothesis: Hypothesis (H) and deg x (P) < d x . 

2. HERMITE REDUCTION 

Let K be a field of characteristic zero, either k or k(x) 
in what follows. Let K(y) be the field of rational functions 
in y over K, and D y be the usual derivation on it. For a ra- 
tional function / G K(y), Hermite reduction [10] computes 
rational functions g and r = a/b in K(y) satisfying 

/ = D v(g) + r, deg y (a) < deg y (b), b is squarefree. (3) 

Horowitz and Ostrogradsky's method [15, 11] computes the 
same decomposition as in (3) by solving a linear system. For 
the details of those methods, see [4, Chapter 2]. 

Lemma 4 If f is proper, a pair (g, r) satisfying (3) for 
proper g,r is unique. 

Proof. This is a consequence of [11, Theorem 2.10] after 
writing rasa sum Y^iLi a i/( x ~ &*) an d integrating. □ 

Lemma 5 Let f be a nonzero rational function in K(y) of 
degree at most n in y, then Hermite reduction on f can be 
performed using 0(n) operations in K . 

PROOF. See [6, Theorem 22.7]. □ 



In contrast, the method of Horowitz and Ostrogradsky takes 
(D(n^) operations in K [6, §22.2]. Thus, Hermite's method 
is quasi-optimal and asymptotically faster than the former. 

From now on, we fix K = k(x) and analyse the complexity 
of Hermite reduction over k(x) in terms of operations in k. 
To this end, we use an evaluation-interpolation approach. 

2.1 Output size estimates 

We derive an upper bound on the bidegrees of g and r 
satisfying (3) by studying the linear system in [11]. 

Analysing Hermite reduction (under (H)) shows the exis- 
tence of A, a G k(x)[y] withdeg y (A) < d~ , deg y (a) < d* and 
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In order to bound the bidegrees of A and a, we reformu- 
late (4) into the equivalent form 



'Dy(A) 
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where Q* D y (Q~)/Q~ is a polynomial in fc[s,y] of bidegree 
at most (d*, d* — 1) by Fact 3. Viewing A and a as polyno- 
mials in fc(a;)[2/] with undetermined coefficients, we form the 
following linear system, equivalent to (5), 



(Hi H2 



= P, 



(6) 



where Hi G fc[ x ]<d* * i H2 G k[x] 1 and A, a, and P 

are the coefficient vectors of A, a, and P with sizes d~ , d*, 
and d y , respectively. Under the constraint of properness of 
A/Q~ and a/Q*, (A, a) is unique by Lemma 4. Then (6) 
has a unique solution, which leads to the following lemma. 

Lemma 6 The matrix (Hi H2) is invertible over k(x). 

As the matrix (Hi H2) is uniquely defined by Q, we call 
it the matrix associated with Q, denoted by H(Q). Let 5 
be its determinant, so that deg x (S) < fi := d x dy + d x d y by 
Fact 2(i). For later use, we also define 8' as the determinant 
oiH{Q* 2 ), so that deg^') < n' := 2d x d* y by Fact 2(i) and 
since (Q* 2 )~ = Q* ■ 

Lemma 7 There exist B, b G k[x,y] with deg y (B) < dy 
and deg y (b) < d y , and such that: 

(ii) deg x (B) < fi - d* + deg x (P) and deg x (b) < 
deg x (P). 

Proof. Applying Cramer's rule to (6) leads to (i). As- 
sertion (11) next follows by determinant expansions. □ 

In what follows, we shall encounter proper rational func- 
tions with denominator Q satisfying Q = Q* 2 . The following 
lemma is an easy corollary of Lemma 7 for such functions. 

Corollary 8 Assuming Q = Q* 2 in addition to Hypothe- 
sis (H), there exist B,b G k[x,y] with deg y (B) and deg y (b) 
less than d* y , and such that 
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Q.2 — \W~Q* J t S'Q* 

(ii) deg x (B) and deg x (b) are bounded by p! — d x + deg x (P) . 



Algorithm HermiteEvallnterp(P, Q) 

Input: P,Q £ k[x,y] satisfying Hypothesis (H). 
Output: (A, a) £ k(x)[y] 2 solving (4). 

1. Compute Q~ := gcd(Q,D y (Q)) and Q* := Q/Q~; 

2. Set A := 2(d x d y + d* y d~) + deg x (P) - mm{d x ,d*}; 

3. Set S to the set of A+l smallest nonnegative integers 
that are lucky for Q; 

4. For each xo £ 5*, compute (Ao,ao) £ k[y] 2 such that 

P(x ,y) = D ( A \ a 
Q(xo,y) v \Q~(x ,y)J Q*{x ,y) 
using Hermite reduction over fe; 

5. Compute (A, a) £ k(x)[y] by rational interpolation 
and return this pair. 



2.2 Algorithm by evaluation and interpolation 

We observe that an asymptotically optimal complexity 
can be achieved by evaluation and interpolation at each step 
of Hermite reduction over k(x). This inspires us to adapt 
Gerhard's modular method [8, 9] to k(x,y). Recall that, by 
Hyp. (H), Q £ k[x,y] is nonzero and primitive over k[x]. 

Definition An element xo £ k is lucky if \c y (Q)(xo) 7^ 
and deg y (gcd(Q(x ,y), D y (Q(x , y)))) = d y . 

Lemma 9 There are at most d x (2d y — 1) unlucky points. 

Proof. Let a £ k[x] be the d~th subresultant w.r.t. y 
of Q and D y (Q). By [9, Corollary 5.5], all unlucky points 
are in the set U — { xo £ k \ a(xo) = 0}. By [9, Corollary 
3.2(H)], deg x {a) < d x {2d* y - 1). □ 

Lemma 10 Let B , b, and 8 be the same as in Lemma 7, and 
let xo £ k be lucky. Then S(xo) 7^ and (B(xo,y),b(xo,y)) 
is the unique pair such that 

P(xo,y) = D ( B ( x o,y) \ + b( x o,y) , 7 s 
Q(xo,y) y \S(x )Q-(x ,y) J S(x )Q*(x ,y)' 

Proof. By the luckiness of xo, deg y (Q(xo,y)) = d y and 
Q(x ,y)~ = Q~(x ,y), so Q(x ,y)* = Q*(x ,y). This im- 
plies rl(Q)(xo,y) — H(Q(xo,y)), which, by Lemma 6, is 
invertible over k(x). Hence 5(xo) 7^ 0, and the evaluation 
at x = xo of the equality in Lemma 7(i) is well-defined. 
Thus, (B(xo,y),b(xo,y)) is a solution of (7). Uniqueness 
follows from Lemma 4. □ 

Theorem 11 Algorithm HermiteEvallnterp in Figure 2 is 
correct and takes 0(d x d y + deg x (P)d y ) ops. 

Proof. Set v to d x (2d* — 1). Lemma 9 implies that the 
A + l lucky points found in Step 3 are all less than A + 
v + 1. By Lemmas 4 and 7(i), A = B/S and a — b/8. By 
Lemma 10, Ao = B(xo,y)/S(xo) and ao = b(xo,y) / &{xo). 
By Lemma 1(H) and since deg a .(5) < [i, it suffices to ratio- 
nally interpolate A and a from values at A + 1 lucky points. 
This shows the correctness. The dominant computation in 
Step 1 is the gcd, which takes 0(d x d y ) ops by [6, Cor. 11.9]. 
For each integer i < \ + testing luckiness amounts to eval- 
uations at xo and computing gcd(Q(xo,y),D y (Q(xo,y))), 
which takes (D(d y ) ops by Fact l(i) and [6, Cor. 11.6]. Then, 
generating S in Step 3 costs (5((A + v + l)d y ) ops. By 
Fact l(i), evaluations in Step 4 take (5((A + l)d y ) ops. For 
each £0 £ •S', the cost of the Hermite reduction in Step 4 is 
0(d y ) ops by Lemma 5. Thus, the total cost of Step 4 is 
<5((A + l)d y ) ops. By Fact l(ii), Step 5 takes £>((A + l)d y ) 
ops. Since A < 2d x d y + deg x (P) and v < 2d x d y , the total 
cost is as announced. □ 

As the generic output size of Hermite reduction is propor- 
tional to Xd y , which is 0((d x d y + deg x (P))d y ), Algorithm 
HermiteEvallnterp has quasi-optimal complexity. 

3. MINIMAL TELESCOPERS 

We analyse two algorithms for constructing minimal tele- 
scopers for bivariate rational functions and their certificates. 



Figure 2: Hermite reduction over k(x) via evaluation 
and interpolation. 

3.1 Hermite reduction approach 

We design a new algorithm, presented in Figure 3, to com- 
pute minimal telescopers for rational functions by basing on 
Hermite reduction. For / = P/Q £ k(x, y) and i £ N, Her- 
mite reduction decomposes D x (/) into 

DUf) = D y ( gi ) + n, (8) 

where gi,fi £ k(x,y) are proper. Since the squarefree part 
of the denominator of D x (f) divides Q* , so does the de- 
nominator of ri. The following lemma shows that (8) re- 
combines into telescopers and certificates; next, Lemma 13 
implies that the first pair obtained in this way by Algorithm 
HermiteTelescoping in Figure 3 yields a minimal telescoper. 

Lemma 12 The rational functions ro, . . . , rd* are linearly 
dependent over k(x). 

Proof. The constraints on ri imply deg y (riQ*) < d* y for 
all i £ N, from which follows the existence of a nontrivial 
linear dependence among the tVs over k(x). □ 

Lemma 13 An integer p is minimal such that X/i=o r l iri 
u for rjo, ■ ■ ■ ,n p £ k(x) not all zero if and only if X^i=o ViDx 
is a minimal telescoper for f with certificate y~]?_ n Jfeflj. 

Proof. Multiplying (8) by 77; before summing yields 

/ p \ p p 
L (f) = D v{y2vi9i) +^2 run for L -^i/,^, 

M = ' i=0 i=0 

where the first two sums are proper. Thus, by Lemma 4, L is 
a telescoper of order p for / with certificate 53i=o 7 l i 9 i ^ anc ^ 
only if 53i=o niri ~ w ith n p 7^ 0. The lemma follows. □ 

3.1.1 Order bounds for minimal telescopers 

Lemmas 12 and 13 combine into an upper bound on the 
order of minimal telescopers for /. 

Corollary 14 Minimal telescopers have order at most d* . 

The bound 6d y is shown in [3] for rational functions of 
the form yD y (Q)/Q with Q £ Ai[a;, 3/]. Apagodu and Zeil- 
berger [2] obtain a similar bound for a class of nonrational 



hyperexponential functions, but their proof does not seem 
to apply to rational functions, as it heavily relies on the 
presence of a nontrivial exponential part. 

We also derive a lower bound on the order of the minimal 
telescoper, to be used as an optimisation at the end of § 3.1.3: 
choosing a lucky xo £ k, next applying Hermite reduction 
in k(y) to D l x (f)(x ,y), yields 



D x (f)(x ,y) = Dy(g 0>i ) +r ,i, 



(9) 



where go t i,ro,i £ k(y) are proper and the denominator of ro,i 
divides Q*(xo,y). Let po be the smallest integer such that 
ro,o, ■ • ■ , ro, PO are linearly dependent over k. 

Lemma 15 A minimal telescoper has order at least po- 

Proof. We first claim that ro,i = ri(xo,y), for n as 
in (8). Note that the squarefree part w.r.t. y of the denomi- 
nator of D x (f) divides Q* for all ieN. By [9, Cor. 5.5], x is 
lucky for the denominator of D x (f) for all i £ N. Then, the 
claim on ro,i follows from Lemma 10 applied to D x (f). Let p 
be the minimal order of a telescoper, then ro, . . . ,r p are lin- 
early dependent over k(x) by Lemma 13. Thus ro,o, • • ■ , ?"o,p 
are linearly dependent over k, which implies po < p. □ 

3.1.2 Degree bounds for minimal telescopers 

To derive degree bounds for gi and n in (8), let S, 8' , p, 
and p! be defined as before Lemma 7, and set p" — p+p' — 1. 

Lemma 16 Let W be in k\x,y] with deg y (W) < d y . Then, 
for all i £ R, there exist B,b £ k[x, y] with both bideg(B) 
and bideg(fe) bounded by (deg^W) + p" , d* — 1), such that 
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Proof. A straightforward calculation leads to 

W \ W 1 WD X (Q*) 



D, 



8 i+2 8' i+1 Q* S^S" 



where bideg(H>) < (deg x (W) + p" , d y - 1). By Corollary 8, 
there exist B, b £ k[x, y] such that 



WD X {Q*) 



Q 



,2 



,i+l 



SB 



81, 



with bideg(B) and bideg(6) bounded by (deg x (W) + p' — 1, 
d* - 1). Setting (B, b) = {—SB, W - 8b) ends the proof. □ 

Lemma 17 For i £ N, there exist Bi,bi £ k[x,y] such that 
B, \ h 



DlU)=D v 



+ 



(10) 



Moreover, bideg(_Bi) < (deg x (P) + p + ip" + (i — l)d*, id* y + 
dy — 1) and bideg(b;) < (deg x (P) + p + i/i" — d^T , d* — 1) . 

Proof. We proceed by induction on i. For i — 0, the 
claim follows from Lemma 7. Assume that i > and that 
the claim holds for the values less than i. For brevity, we 
set 7 = deg x (P) + p, Fi-i = B 1 . 1 /(S l 8"~ 1 Q"- 1 Q-), and 
Gi-i = bi-i/(8 % S n ^ 1 Q*). The induction hypothesis implies 

Dl(f) = D y D x (F l - 1 ) + D x (Gi-i), 

with bidegree bounds on Bi-i and 6»_i. Fact 3(i) implies 
that Q :— Q* D X (Q~)/Q~ is in k[x,y], with bideg(Q) < 



Algorithm H e r m i teTe I esco p i n g ( / ) 

Input: / = P/Q e k(x,y) satisfying Hypothesis (H). 
Output: A minimal telescoper L £ k[x](D x ) with cer- 
tificate g £ k(x,y). 

1. Apply HermiteEvallnterp to / to get (go,ao) such 
that / = Dy(go) + a /Q* . If a = 0, return (l,flo)- 

2. For i from 1 to deg y (Q*) do 

(a) Apply HermiteEvallnterp to ~ai~iD x (Q*)/Q* 2 
to express it as D y (§i) + cii/Q* . 

(b) Set gi = D x (g i -i)+g 1 and a, = D x {a i - 1 ) + a t . 

(c) Solve X!j=o r ?j a J = f° r % e k(x) using [19]. 
If there exists a nontrivial solution, then set 
{L,g) := ,.'//':•>], :;'/..'/,)• and break. 

3. Compute the content c of L and return (c~ 1 L, c~ 1 g). 



Figure 3: Creative telescoping by Hermite reduction 

(d* — l,d*). Hence D x (l/Q~) = —Q/Q. This observation 
and an easy calculation imply that 

D (F ) = Bj-i 

where _Bi_i £ fc[i,j/] and deg ;c (Bi_i) < deg^Bi-i) + p" + 
d*. Furthermore, by Lemma 16 there are Bi,bi £ k[x,y] 
with bidegrees at most (deg a .(6i_i) + p" ,d* — 1), such that 



8 l + 1 8' l Q* ) 

Setting Bi = + BiQ* 1 ~ 1 Q~ and 6i = 6i, we arrive 

at (10). It remains to verify the degree bounds. The induc- 
tion hypothesis implies that both deg x (Bi) and deg x (bi) are 
bounded by j + ip" — d x . It follows that deg x (BiQ* 1 ~ 1 Q~) 
is bounded by 7 + ip" + (i — l)d*. Similarly, deg x (Bi-i) is 
bounded by 7 + ip" + (i — l)dj, and so is deg a .(Bi). The 
bounds on degrees in y are obvious. □ 

We next derive degree bounds for the minimal telescop- 
ers obtained at an intermediate stage of HermiteTelescoping; 
refined bounds on the output will be given by Theorem 25. 

Lemma 18 Under (H'), Step 2(c) of Algorithm Hermite- 
Telescoping computes a minimal telescoper L £ k[x]{D x ) 
with order p and a certificate g £ k(x,y) for P/Q with 
deg x (L) £ 0(d x d y p 2 ) and bideg(g) £ 0(d x d y p 2 ) x 0(d y p). 

Proof. By Lemma 13, we exhibit a minimal telescoper 
by considering the first nontrivial linear dependence among 
the ai's in (10). Let M be the coefficient matrix of the 
system in (r)i) obtained from X^i > =o 7 ? i<Ii = 0- By Lemma 17, 
M is of size at most (p + 1) x dy and with coefficients of 
degree at most a :— d x + p + pp" — d x in x. Hence, there 
exists a solution (770, . . . , rjp) £ fc[a;] p+1 of degree at most op 
in x by Fact 2(ii). Since p, p" £ 0(d x d y ) and d* < d y , the 
degree estimates of L and g are as announced. □ 

3.1.3 Complexity estimates 

We proceed to analyse the complexity of the algorithm in 
Figure 3 and of an optimisation. 



Theorem 19 Under Hyp. (H'), Algorithm HermiteTelescop- 

ing in Figure 3 is correct and takes 0(p^ +1 d x d 2 ) ops, where 
p is the order of the minimal telescoper. 

Proof. The formulas in Step 2(a) create the loop invari- 
ant D x (f) = Dy(gi) + ca/Q* . Correctness then follows from 
Lemmas 12 and 20. Step 1 takes (D(d x d 2 ) ops by Theo- 
rem 11 under (H'). By Lemma 17, deg a .(— Oi-iD x (Q*)) £ 
0(id x d y ). So the cost for performing Hermite reduction on 
-a i - 1 D x (Q*)/Q* 2 in Step 2(a) is d(id x d 2 y ) ops by Theo- 
rem 11. The bidegrees of <?; and ai in Step 2(b) are in 
0(id x dy) x 0(id y ) by Lemma 17. Since adding and differ- 
entiating have linear complexity, Step 2(b) takes (D(i 2 d x d 2 ) 
ops. For each i, the coefficient matrix of Y^=o r l3 a j ~ 
in Step 2(c) is of size at most (i + 1) X d y and with coef- 
ficients of degree at most deg x (ai) £ 0{id x d y ). Moreover, 
the rank of this matrix is either i or i + 1. Then, Step 2(c) 
takes 0(i" d x d y ) ops by Fact 2(iii). Computing the content 
and divisions in Step 3 has complexity 0(d x d y p 3 ). If the 
algorithm returns when i = p, then the total cost is in 
p p 

Y otfdxdl) + Y o{rd x d 2 y ) c 6{ P " +1 d x dl) ops, (ii) 

which is as announced. □ 

An optimisation, based on Lemma 15, consists in guessing 
the order p so as to perform Step 2(c) a few times only: As 
a preprocessing step, choose xo £ k lucky for Q, then detect 
linear dependence of {ro,o, • ■ • , ro,j} in (9). The minimal j 
for dependence is a lower bound po on p. So Step 2(c) is 
then performed only when i > po. In practice, the lower 
bound po computed in this way almost always coincides with 
the actual order p. So normalising the g^s becomes the 
dominant step, as observed in experiments. We analyse this 
optimisation by first estimating the cost for computing po. 

Lemma 20 Under Hypothesis (H'), computing a lower or- 
der bound po for minimal telescopers takes O(d x d y po) ops. 

Proof. Since differentiating has linear complexity, the 
derivative D' x (f) takes 0(i 2 d x d y ) ops. By Fact l(i), the 
evaluation D x (f)(xo,y) takes as much. The cost of Hermite 
reduction on D x (f)(xo,y) is (D(id y ) ops by Lemma 5. By 
Fact 2 (Hi) with d = 1, computing the rank of the coefficient 
matrix of YTj=o r lj r o^> with r j as in (9), takes 0(d H i" _1 ) 
ops. Thus, the total cost for computing a lower bound on po 
is Ei=o 0{i 2 d x d y ) e 6{d x d y pl) ops. □ 

Corollary 21 For runs such that po = p — 0(1), the previ- 
ous optimisation of HermiteTelescoping takes <D{p' i d x d 2 l ) ops. 
Proof. In view of Lemma 20, the estimate (11) becomes 

d{d x d y pl) + 2Xo d(i 2 d x d 2 y ) + Ef =po 6(i"d x d 2 y ), which is 

0(p 3 d X dy) + 0((P - P 0)p U 'd X d 2 y) ops, whence the result. □ 

3.2 Almkvist and Zeilberger's approach 

We analyse the complexity of Almkvist and Zeilberger's 
algorithm [1] when restricted to bivariate rational functions. 
In order to get a telescoper whose order p is minimal, the re- 
sulting algorithm, denoted RatAZ, solves (1) for increasing, 
prescribed values of p until it gets a solution (rjo, ■ ■ ■ , r\ p , g) £ 
k(x) p+1 x k(x, y) with the r/i's not all zero. For the analysis, 
we start by studying the parameterisation of the differential 
Gosper algorithm of [1] under the same restriction to k(x, y). 



Definition ([9]) Let K be a field and a, b £ K[y] be non- 
zero polynomials. A triple (p,q,r) 6 -Kly] 3 is said to be a 
differential Gosper form of the rational function a/b if 

- = £l>M. + 1 and gcd(r, q - rDJr)) = 1 for all t6N. 
b p r 

For hyperexponential /, a key step in [1] is to compute a 
differential Gosper form of the logarithmic derivative of F = 
X/i=o T H^xif)t w h ere the iji's are undetermined from k(x). 
In the analogue RatAZ, this form is predicted by Lemma 22 
below, which is a technical generalisation of a result by 
Le [12] on F when / has a squarefree denominator. 

Write Q — t(y)T(x,y), splitting content and primitive 
part w.r.t. x. By an easy induction, D x (f) = Ni/(QT* t ) 
for Ni e k[x,y]. For this section, set F — 5Z£ =0 »?i-Di(/), 
N = E£=o ViNiT* p -\ and H = -D y (Q)/Q~ -pt*D y (T*). 

Lemma 22 If F is nonzero, the triple (N, H,Q*) is a dif- 
ferential Gosper form of D y (F) /F . 

Proof. First, observe F = N/(QT* P ) and Q* = t'T*. 
Next, D y (F)/F = D y (N)/N - D y (Q)/Q - P D y (T*)/T* is 
D y (N)/N + H/Q*. There remains to prove gcd(Q*,H - 
tD v (Q*)) = 1, for any r £ N. Recall that the squarefree 
part Q* of Q is the product Q1Q2 ■ ■ ■ Qm and that Qi de- 
notes Q* /Qi. By Fact S(ii), 

m 

Z:=H- rD y (Q*) = -pt*D y (T*) - ]T(i + r)Q t D y (Q t ). 

i=l 

If Qj divides t* , Z reduces to —(j + r)QjD y (Qj) modulo Qj. 
If not, it reduces to -{3 J rT)QjD y {Q j )-pt*[D y {Q j )T*/Q j ), 
which rewrites to —(j+T+p)QjD y (Qj) modulo Qj. In both 
cases, Z is coprime with Q* , as j > 0, r > 0, and p > 0. □ 

By another induction, we observe bideg(Ai) < (deg a .(P) + 
ideg x (T*) - i,d y + ideg y (T*) - 1), so that bideg(Ar) < 
(degJP) +pdeg m (T*) - p, d y +pdeg y (T*) - 1). 

The next step in RatAZ is, for fixed p, to reduce (1) by 
the change of unknown g — z/(Q~T* p ), so as to determine 
all (r?i) £ k(x) p+1 for which the differential equation in z 

p 

Y^ViNiT*"^ = Q'D y {z) + (D y (Q*) + H) z (12) 

i=0 

has a polynomial solution in k(x)[y]. For later use, we recall 
the following consequence of [9, Corollary 9.6]. 

Lemma 23 Let a,b £ K[y] be such that /3 — — lc y (b) / lc y (a) 
is a nonnegative integer and deg y (b) = deg y (a) — 1. Let 
c £ K[y] be such that (3 > deg H (c) — deg H (a) + 1. If u is a 
polynomial solution of aD y (z) + bz = c, then deg y (u) < /3. 

The following lemma generalises [12, Lemma 2] to present a 
degree bound for z. 

Lemma 24 If u £ k(x)[y] is a solution of (12) for (qi) £ 
k{x) p+ , then deg y (u) is bounded by /3 = dy + pdeg y (T*). 

PROOF. Let a = Q* and 6 = D y (Q*) + H. By the defi- 
nition of H, b = -Q*D v (Q-)/Q- - pt*D y {T*) '. Fact 3(i) 
implies that lc^(6) = — (d y + pdeg y (T*)) \c y (a). Therefore, 
= -lc y (b)/lc y (a) = d y +pdeg y (T*). As deg H (A0 < 
d y + pdeg y (T*) and d y = d* y + d~ , /3 > deg y (N) - d y + 1. 
The lemma holds by Lemma 23. □ 



Algorithm RatAZ(/) 

Input: / = P/Q e k(x,y) satisfying Hypothesis (H). 
Output: A minimal telescoper L £ k[x](D x ) with cer- 
tificate g £ k(x,y). 

1. Compute Q~ = gcd(Q, Q* = Q/Q~, and 
T, T* primitive parts of Q, Q* w.r.t. x, respectively; 

2. Set (N,N,/3,H) to (P,P,d y ,-Q*D y (Q)/Q); 

3. For £ = 0, 1, . . . do 

(a) Set z to ^2j =Q Zjy : ', extract the linear system 

Ai (rji Zj~) T = from (12) (for p — £) and com- 
pute a basis S of the null space of Ai by [19]. 

(b) If S contains a solution (770, . . . , r\i , s) such that 
rjo, . . ■ , rji are not all nonzero, then set (L, g) := 
{IlLo ViDl,s/(Q-T* e )), and go to Step 4; 

(c) Update N := D X (N)T* - N(T*D x (T)/T + 
iD x {T*)), N := NT* + m+1 N, f3 ~ /3 + 
deg H (T*), and H := H — t*D y {T*). 

4. Compute the content c of L and return (c _1 L, c~ 1 g). 



Figure 4: Improved Almkvist— Zeilberger algorithm 

We end the present section using the approach of Almkvist 
and Zeilberger to provide tight degree bounds on the outputs 
from Algorithms HermiteTelescoping and RatAZ. 

Theorem 25 Under Hypothesis (H'), there exists a mini- 
mal telescoper L £ k[x]{D x ) with certificate g £ k(x,y) with 
deg x (L) € 0(d x d y d y ) andbideg(g) € 0{d x dyd* y )xO{d y d* v ). 

Proof. By Corollary 14, there exists a smallest p £ N 
at most d v , for which (1) has a solution with the n's not 
all zero. For this p, we estimate the size of the polynomial 
matrix Ai derived from (12) by undetermined coefficients. 
By the remark on iV after Lemma 22, we have bideg(iV) < 
(n x , n y ) where n x := d x +pdeg x (T*) — p £ G(pd x ) and n y := 
d y + pdeg y (T*) — 1 £ O(pdy). The matrix Ai contains two 

blocks Mi £ fcW^ x +1)X(P+1> an d M2 £ k[x] { ^+ 1)x(0+1) , 
where /3 £ 0(pd y ) is the same as in Lemma 24. By the 
minimality of p, the dimension of the null space of Ai is 1. 
So there exists u £ fc[a;] n! ' + 1 with coefficients of degree at 
most n x {p + 1) + d x (/3 + 1) £ 0(d x d y d v ) in x such that 

Ai (rj 2) T = 0, which implies degree bounds in x for L and g. 
The degree bound in y for g is obvious. □ 

We now analyse the complexity of the algorithm in Fig. 4. 

Theorem 26 Under Hypothesis (H'), Algorithm RatAZ in 
Figure 4 is correct and takes <D(d x dy p" +2 ) ops, where p is 
the order of the minimal telescoper. 

Proof. By the existence of a telescoper, Corollary 14, 
and Lemma 24, the algorithm always terminates and returns 
a minimal telescoper L, of order p at most dy. Gcd computa- 
tions dominate the cost of Steps 1 and 2, which take 0(d x dy) 
ops. For each £ £ N, the dominating cost in Step 3 is com- 
puting the null space of Ai. Let n y — dy +£deg y (T*) — 1 £ 
0{£d y ) and n x = d x + £deg x {T") £ 0(£d x ). By the same 
argument as in the proof of Theorem 25, the matrix M 



is of size at most (n y + 1) x (£ + (3 + 2) and with coeffi- 
cients of degree at most n x . Let r be the rank of Ai, which 
is either £ + f3 + 2or£ + f3 + l by construction. Thus, 
a basis of the null space of Ai can be computed within 
0{n x (n y + 1)(£ + 13 + 2)r"^ 2 ) ops by Fact 2(m). Since 
13 £ 0(£d y ), 6{n x (n y + 1)(£ + (3 + 2)r"" 2 ) is included in 
0(d x dy£" +1 ). Since Step 3 terminates at £ = p, the total 
cost of the algorithm is ~Y^ P l=0 d x dyt }+1 ops. This is within 
the announced complexity, 0(d x d y p" +2 ) ops. □ 

Corollary 27 Algorithms HermiteTelescoping and RatAZ in 

Fig. 3 and 4 both output the primitive minimal telescoper L 
together with its certificate g, which satisfy deg D (L) < d* y , 
deg J .(L),deg a .(s) £ 0(d x d y d y ), and deg y (g) £ 0(d y d*). 

Proof. Both algorithms output the primitive minimal 
telescoper, as they compute a minimal telescoper at an inter- 
mediate step, and owing to their last step of content removal. 
Bounds follow from Corollary 14 and Theorem 25. □ 

4. NONMINIMAL TELESCOPERS 

Here, we discard Hypothesis (H) and trade the minimality 
of telescopers for smaller total output sizes. To this end, we 
adapt and slightly extend the arguments in [13] and [3, §3]. 

Given / = P/Q £ k(x,y) of bidegree (d x ,d y ), our goal 
is to find a (possibly nonminimal) telescoper for /. It is 
sufficient to find a nonzero differential operator A(x, D x ,D y ) 
that annihilates /. Indeed, any A £ k[x]{D x , D y ) \ {0} such 
that A(f) - can be written A = D r y {L + D y R), where L 
is nonzero in k[x](D x ) and R £ k[x]{D x , D y ). If r = 0, then 
clearly L is a telescoper for /; otherwise, A(f) — yields 
Hf) = D v (-R(f) - Efco^2/ i+1 ) for some a t £ k(x), 
which implies that L is again a telescoper for /. Moreover, in 
both cases, deg x (L) < deg x (A) and deg Cx (L) < deg Dx (A). 
Furthermore, for any (i, j, £) £ N 3 , a direct calculation yields 

^DiDt(f) = j^j. (13) 

where Hi, j:t £ k[x, y] and deg^. (H ii:j , e ) < (j + £ + l)d x + i - j 
and deg y (Hij t e) < (j + £+ l)d y — £. From these inequali- 
ties, we derive the size and complexity estimates in Figure 1 
(bottom half), using two different filtrations of k[x]{D x , D y ). 

Lipshitz's filtration ([13]). Let F v be the k- vector space 
of dimension f„ := ("^ 3 ) spanned by {x l D x D y \ i+j + l < 
v}. By (13), Fv(f) is contained in the vector space of di- 
mension g„ := ((1/ + l)d x + v + l) {{v + l)d y + 1) spanned 

by { I i < {v + l)d x + v, j < (v + l)d y }. Choos- 

ing v — 6(d x + l)(d y + 1) yields f„ > g„; therefore, there 
exists A in k(x, D x , D y ) \ {0} with total degree at most 
Q(d x + V)(d y + 1) in x, D x , and D y that annihilates /. More- 
over, A is found by linear algebra in dimension 0((d x d y ) 3 ). 
A better filtration ([3]). Instead of taking total de- 
gree, set F K<V to the fc-vector space of dimension f Kj „ := 
(ft + generated by {VDjD* | i < k, j + i < v}. 

By (13), F KtW (f) is contained in the vector space of dimen- 
sion g Ki „ := ((// + l)d x + n + l)((v + l)d y + 1) spanned by 

{ ^rr I i < {v + l)d x + k, j < {v + l)d y ). Choosing 
K = 3d x d y and v — 6d y results in f Kj „ > g K .„. This implies 
the existence of A in k{x, D x , D y ) \ {0} with total degree at 
most 6d y in D x and D y and degree at most 3d x d y in x that 
annihilates /. Again, A is found by linear algebra over k, 
but in smaller dimension 0(d x d y ). 
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Table 1: Creative telescoping on random instances 
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Timings in ms for algorithms in Table 3 (stopped after 30 min). 

5. IMPLEMENTATION AND TIMINGS 

We implemented in Maple 13 all the algorithms described; 
as we used Maple's generic solver SolveTools : -Linear, all 
of our implementations are deterministic. 

The evaluation-interpolation algorithm HermiteEvallnterp 
for Hermite reduction (Fig. 2) does not perform well, mainly 
because Maple's rational interpolation routines are far too 
slow. We thus implemented Algorithm HermiteReduce (orig- 
inal version) in [4, § 2.2] (carefully avoiding redundant ex- 
tended gcd calculations), and noted that it performs better. 

We then implemented a variant of Algorithm HermiteTe- 
lescoping in Figure 3, using HermiteReduce in place of Her- 
miteEvallnterp, and including the optimisation at the end of 
§3.1.3, refined by additional modular calculations. 

For a rational function, Algorithm HermiteTelescoping re- 
turns the minimal telescoper L and the certificate g. The 
algorithm separates the computation for L from that for g. 
Indeed, g is formed by the coefficients of L, go, the §i and 
their derivatives given in Figure 3. This feature enables us 
to either return the certificate g as a sum of unnormalised 
rational functions, or a normalised rational function. 

A selection of timings by this implementation and oth- 
ers are given in Table 1; our code, the full table, as well as 
the random inputs are given in [22]. For our experiments, 
we exhaustively considered all 49 bidegree patterns in fac- 
torisations of denominators Q% ■ ■ ■ Q™ (m < 5) that add up 
to bidegree (5,5), and generated corresponding random de- 
nominators, imposing the integers of the expanded forms to 
have around 26 digits. Numerators were generated as ran- 
dom bidegree-(5,5) polynomials with coefficients of 26 digits. 

Application to diagonals. The diagonal of a formal power 
3 k[[x, y]] is defined to be the 
For a D-finite power se- 
ries /, it is known to be D-finite [13], and it is even algebraic 
for a bivariate rational function / £ k(x,y) n fc[[x,y]] [18, 
§6.3]. A linear differential operator L £ k(x){D x ) that an- 
nihilates A(/) can then be computed via rational-function 
telescoping, owing to the following classical lemma from [13] . 

Lemma 28 Any telescoper for f(y, ^)/y annihilates A(/). 

By this lemma, it suffices to compute a telescoper without its 
certificate to get an annihilator. Algorithm HermiteTelescop- 
ing is suitable for this task, since it separates computation of 
telescopers and certificates. Alternatively, for / = P/Q, we 
can compute an annihilator of A(/) either as the differen- 
tial resolvent of the resultant Res H (Q, P — rD y Q), or simply 
guess it from the first terms of the series expansion of A(/). 

We compare the various algorithms on an example bor- 
rowed from [5] (timings of execution are given in Table 2): 



series / = E^>o fi,i x 'V m k i\ x 
power series A(/) := Y^Zo fi.i x 



f 



1 



1 — x — y — xy(l — x d ) 



where d £ N. 



(14) 



Table 2: Computation of the diagonals of (14) 
Timings in ms by creative telescoping of f(y,x/y)/y (upper half) 
or f(y/x,x)/x (second half). Algorithms listed in Table 3. 

AZ DETools [Zeilberger] 
Abr AZ with Abramov's denominator bound by option gosper_f ree 
RAZ Algorithm RatAZ of Fig. 4, with lower-bound prediction 

HI our Hermite-based approach, without certificate normalisation 

H2 HI, but with normalised certificate 

HQ RAZ, solving (1) by Horowitz— Ostrogradsky 

EI HI with evaluation and interpolation for calculations over k(x) 
MG Mgfun's creative telescoping for general D-finite functions 
RR telescoper computation by resultant and differential resolvent 
GHP telescoper guessing by diagonal expansion and Hermite— Padc 

Table 3: List of the algorithms for the experiments 

6 



All computer calculations have been performed on a Quad- 
Core Intel Xeon X5482 processor at 3.20GHz, with 3GB of 
RAM, using up to 6.5GB of memory allocated by Maple. 
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